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I.  INTRODUCTION 


Since  Che  formulation  of  Che  governing  equeciona  and  pioneering  work  by 
researchers  on  solid  rocket  nozzle  flows  in  Che  1960s  (e.g..  Refs.  1  and  2), 
Che  CheoreCical  approach  Co  gas-particle  flows  has  received  increasing 
attention  by  the  workers  in  the  propulsion  industry  who  are  involved  with  the 
development  of  solid  rocket  motors  (SRMs).  The  knowledge  gained  from 
theoretical  two-phase  SRM  flowfield  studies  will  lead  to  better  predictions  of 
motor  performance,  range,  and  payload  delivery,  improvement  of  motor  in-depth 
thermal-structural  analysis,  and  reduction  in  motor  design,  development,  and 
test  risks  and  cost.  A  complete  description  of  the  SRM  flow  field  necessi¬ 
tates  Che  consideration  of  combustion  chamber,  igniter,  throat,  and  exhaust 
exit  cone  internal  fully  coupled  chemically  reacting  viscous  two-phase  flow 
solution,  which  is  by  no  means  an  easy  task.  Only  the  analysis  of  flow  under 
idealized  conditions  is  possible  at  the  present  time.  In  the  subsonic-tran¬ 
sonic  regime,  Ref.  3  provides  a  solution  method  for  a  two-phase  inviscid  flow 

inside  Che  chamber  and  nozzle  of  arbitrary  configurations.  For  the  supersonic 

4 

nozzle  and  exit  cone  flow,  the  widely  used  computer  program  in  the  industry 
utilizes  the  method  of  characteristics.  There  is  nothing  wrong  with  applying 
the  method  of  characteristics  to  the  two-phase  supersonic  flow  study,  except 
that  the  extension  of  the  method  Co  three-dimensional  space  is  not  straight¬ 
forward.  The  inability  Co  resolve  flowfields  containing  shock  waves  is 
another  drawback  of  adopting  the  method  of  characteristics,  since  almost  any 
disturbance  in  a  supersonic  stream  will  produce  one  or  more  shock  waves.  A 
finite  difference  approach  has  recently  been  applied  to  nozzle  and  exhaust 
plume  flows  including  gas-particle  interactions.^ 

All  the  works  cited  above  have  been  restricted  to  the  flow  in  one-  or 
two-dimensional  space.  Confronted  with  asymmetric  nozzle  mechanizations,  such 
as  the  canted  Titan  SRM,  fluid-bearing  thrust  vector  control  (TVC)  inertial 
upper  stage  (IUS)  motor,  and  the  flexible  bearing  TVC  space  shuttle  SRM,  the 
rocket  nozzle  designer  often  appears  in  an  unfavorable  position.  Currently, 
there  is  no  analytical  tool  available  in  this  regard.  The  scarce  data 
obtained  from  expensive  test  measurements  is  usually  the  only  source  of  design 
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basis.  It  is,  therefore,  highly  desirable  to  have  the  theoretical  capability 
analysing  a  two-phase  flow  inside  a  three-dimensional  nozzle. 

In  this  report,  the  three-dimensional  one-phase  formulation  of  Ref.  6  is 
extended  to  include  the  momentum  and  energy  transfer  between  gas  and  particle 
phases  in  supersonic  nozzle  flows.  The  condition  of  constant  entropy  on  the 
boundary  surface  is  inadequate  for  the  two-phase  flow,  even  in  the  region 
where  a  shock  does  not  interact  with  the  boundary.  Moreover,  the  expansion 
and  compression  wave  structures  in  a  two-phase  nozzle  flow  are  observed  to 
differ  from  those  in  a  clean  gas  nozzle  flow.  The  effect  of  different 
particle  size  and  gas  specific  heat  ratio  on  the  overall  two-phase  flow 
behavior  is  investigated  through  the  calculation  of  a  Mach  3  inlet  flow  inside 
a  rounded  square  three-dimensional  supersonic  nozzle.  Additional  two-phase 
flow  solutions  using  similar  inlet  flow  conditions  for  various  nozzle 
configurations  are  also  presented.  Isometric  projection  of  three-dimensional, 
contour  plots  is  used  for  concise  interpretation  of  the  computed  results  and 
for  visualization  of  flow  structures  inside  three-dimensional  nozzles. 


II.  FORMULATION 


Normalized  by  Che  gas-phase  stagnation  state  corresponding  to  the 
condition  at  the  inlet  plane,  the  governing  equations  written  in  weak 
conservative  form  for  a  steady  three-dimensional  two-phase  flow  take  the 
following  form: 
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with  friction  term 


A. 

J 


1  ill  _i l 

2  m.T 2.  V 

j  j  maxi 


(2  J 


and  energy  exchange  term 


»j  -  2r  iJj  •  ^  -  sc<yi>  -  gc  «.<tJ  -  i4)] 
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where 


N  76  f.P  , 
uj  j  r’ 
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The  momentum  transfer  parameter  f;  is  defined  as 


fi  ■  VS 


Stokes 


(4) 


where  CQ  is  the  particle  drag  coefficient  given  in  Ref.  7  and  CQ 


Stokes 


24/Re ^ 


The  heat  transfer  parameter,  particle  Nusselt  number,  is  taken  as 

N  .  -  2  +  0.459  R°:55  P0'33 
uj  ej  r 

The  particle  Reynolds  number  is  based  on  the  relative  speed 
IAq.1  *  J  (u-u,)2  *  (v-v.)2  +  (w-w.)2 

J  J  J  J 

and  is  defined  as  follows: 
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The  gas  viscosity  is  evaluated  from 


“ti  (I/Tu> 


(7) 


As  shown  in  Fig.  1,  the  geometry  of  the  three-dimensional  nozzle 
cross-sectional  profiles  is  assumed  to  have  at  least  one  plane  of  symmetry, 
and  the  nozzle  wall  radial  coordinate  is  described  by  ry(0,  z).  The 
transformation  of  the  physical  irregular  region  of  interest  into  a  portion  of 
a  unit  circle  with  the  clustering  of  grid  points  in  the  region  of  greatest 
slope  change  is  carried  out  in  the  same  fashion  as  that  in  the  one-phase 
flow.6  The  transformation  relationship  is 


c  - 
5  “ 


2  -1 

—  tan  Oc^tan  [(a  +  bn)])7T  +  a^ 


A 


where  H  *  r/ry  *n<*  constants  k^,  a,  b,  a^,  and  b^  are  parameters 
used  for  the  grid  clustering  control  such  that  for  0  <  ^  <  1 


*  *  “1/2,  b  ■  1/2,  a1  «  1,  bx  -  1  grid  points  are 
clustered  near  n  ■  0. 


a  ■  0,  b  ■  1/2,  a^  ■  0,  b^  ■  1  grid  points  are 
clustered  near  n  ■  1. 

*  “  “1/2,  b  ■  1,  a^  ■  1,  b^  ■  2  grid  points  are 
clustered  near  n  ■  0  and  n  ■  1. 

The  cross-sectional  region  of  computation  for  the  cylindrical  coordinates 
in  Ref.  6  was  restricted  to  90°;  whereas  in  this  two-phase  study,  the 
computational  region  can  be  90°  or  expanded  to  180°.  In  the  transformed 
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Fig.  1.  Three-Dimensional  Nozzle  Coordinate  System 


space  designated  by  a  (z,£,£)  coordinate  system,  the  governing  equations 
take  the  following  form: 


where  the  vectors  E, 
follows: 


3E  3F  3G  . 

3z  +  3;  +  3^  "  H  “  0 


F,  G,  and  H  are  related  to  E,  F,  G, 


(9) 

and  H  of  Eq.  (1)  as 


E  -  i 

F  “  Cn  1*^  *  ^  *  Er^l 
G  «  G 

H  -  H  ~  [CI^  ♦  £nznl  -  Czn  [Fr^  ♦  Gr|g  ♦  £nj  (10) 

The  subscripted  variables  indicate  partial  derivatives. ^ 


III.  SOLUTION  METHOD 


The  weak  conservative  formulation,  Eq.  vC),  is  a  hyperbolic  type  for 

supersonic  flow  along  the  main  flow  z  direction  and  is  solved  by  the 

q 

MacCormack  finite  difference  scheme.  The  calculational  procedure  which 

utilizes  the  stability  analysis  described  in  Ref.  6  and  a  tangency  condition 

for  both  gas  and  particle  phase  at  the  boundary  surface  are  applied  to  the 

two-phase  flow  problem.  Since  two-phase  flow  is  neither  homentropic  in  the 

field  nor  isentropic  at  the  boundary  with  or  without  the  presence  of  shock 

interactions,  the  Abbett  scheme^  for  determining  gas  phase  boundary  flow 

variables  is  modified  to  account  for  entropy  change  at  the  boundary.  The  flow 

9 

variables  obtained  from  the  application  of  predictor  and  corrector  steps  of 
the  finite  difference  method,  in  general,  will  not  satisfy  the  tangent  flow 
condition  on  the  nozzle  surface  in  an  inviscid  analysis,  and  the  Abbett  scheme 
amounts  to  correcting  the  pressure  through  the  use  of  a  simple  wave  to  rotate 
the  flow  vector  to  the  surface  tangent  direction. 

In  internal  flows,  the  small  correction  angle  on  the  nozzle  wall  has  an 
opposite  sign  from  that  of  external  flows  in  Ref.  11.  Instead  of  applying  a 
constant  entropy  condition  to  evaluate  the  gas  density  at  the  boundary,  in 
this  study  Che  finite  boundary  entropy  change  obtained  from  that  of  the 
adjacent  field  points  is  utilized  for  Che  evaluation  of  gas  density  at  the 
boundary.  The  similar  concept  of  using  the  adjacent  field  point  to 
approximate  the  wall  entropy  change  in  a  shock-capturing  approach  was 
introduced  in  Ref.  12  for  single-phase  flows  and  extended  to  multiphase  flows 
in  Ref.  13.  The  three  velocity  components  of  the  gas  phase  are  computed  from 
a  surface  tangency  condition  similar  to  that  given  in  Refs.  6  and  11.  For 
solid  phase,  the  particle  temperature  and  density  are  decoded  from  the 
conservative  variables,  which  are  evaluated  from  a  one-sided  difference 
scheme,  and  the  surface  tangency  condition  is  again  applied  to  three  particle 
velocity  components.  In  reality,  particles  would  either  "stick"  to  or 
"reflect"  from  the  wall^  in  a  complicated  manner,  which  accounts  for  nozzle 
wall  erosion  and  is  still  not  well  understood. 
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An  idealized  situation  where  the  normal  component  o£  the  particle  velocity 
is  zero  has  been  utilized  in  this  study.  A  similar  calculational  procedure 
which  evaluated  the  flow  variables  at  singular  centerline  from  an  averaging 
process  shown  in  Ref.  6  is  applied  to  this  two-phase  study. 


IV.  COMPARISON  WITH  SPP  RESULTS — AXISYMMETRIC  NOZZLE 


To  demonstrate  the  advantage  of  the  finite  difference  approach  to 

two-phase  supersonic  nozzle  flows  and  to  establish  the  credibility  of  the 

CY3D2P  (three-  dimensional  two-phase  flows  in  cylindrical  coordinates)  code, 

developed  under  this  study,  the  results  of  the  calculation  for  an  axisymmetric 

4 

IUS  nozzle  are  compared  with  those  of  the  well-known  SPP  program.  Figure  2 
shows  the  nozzle  and  exit  cone  geometry  for  an  IUS  small  motor.  The  gas  and 
particle  properties  are  the  same  as  those  given  in  Ref.  3  as  follows: 


C  -  1.88  KJ/kg-K  (0.45  Btu/lb  -°R) 
p  m 

C.  *  1.34  KJ/kg-K  (0.32  Btu/lb  -°R) 

J  « 

H  ,  -  8. 444x10 ~5Pa.s  (5.67xl0~5lb  /ft-s) 
tl  m 

5.  -  3203.69  kg/m3  (200  lb  /ft3) 

J  ® 

Pr  -  0.269  0  *  30%  A  -  0.65 

,  r .  -  2.5u  Y  *  1.19 

J 

The  flow  variables  at  Z/rt  “  2.32  obtained  from  the  SPP  code  (1981) 
and  shown  in  Fig.  3  are  utilized  as  initial  data  for  CY3D2P  calculation. 
Figure  4  compares  the  results  from  the  two  programs  at  Z/ r£  ”  10.  The  same 
number  of  initial  data  points  (■  30)  is  used  for  both  calculations.  In 
general,  fairly  good  agreement  between  the  results  of  the  present  finite 
difference  approach  and  those  of  the  updated  1981  SPP  code  is  obtained.  The 
SPP  code  of  the  1975  version  contained  inadequate  centerline  calculation 
procedures  and  would  result  in  erroneous  flow  variables  at  the  singular 
centerline.  For  detailed  flowfield  comparison,  the  updated  1981  SPP  code 
should  be  used. 


'  •*+  •  -  ♦  v'f  k .  ,>y  -  • . 
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^8*  IUS  Small  Motor  Exit  Cone  Configuration 


For  Che  nozzle  geometry  and  gas-particle  properties  considered,  che 
particle  phase  in  the  SPP  code  is  still  undergoing  a  liquid-to-solid  phase 
change  state  and  maintains  at  a  constant  particle  solidification  temperature. 
Mo  provision  has  been  made  in  the  CY3D2P  program  for  particle  phase  change, 
and  a  lower  particle  temperature  associated  with  flow  expansion  in  the  exit 
cone  is  observed  in  the  finite  difference  spproach  than  that  of  SPP  code. 

Since  different  particle  drag  coefficient  formulas  are  used  in  the  two  methods, 
no  such  good  agreement  between  the  results  from  the  method  of  characteristics 
and  those  from  the  present  finite  difference  technique  in  the  two-phase  flow 
calculations  is  expected  as  chat  in  the  one-phase  flow  given  in  Ref.  6.  More¬ 
over,  the  present  results  show  some  "wiggles"  near  the  region  of  particle  free 
zone  at  2/r ^  *  10,  which  is  a  consequence  of  adopting  abrupt  change  of 
particle  flow  variables  for  the  distinct  limiting  particle  streamline  at  the 
initial  data  plane,  the  axial  flow  Mach  numbers  for  different  particle  sizes, 
illustrated  in  Fig.  5,  show  that  increasing  particle  size  does  not  necessarily 
result  in  two-phase  flow  approaching  one-phase  flow  solution  in  a  supersonic 
nozzle.  An  obvious  way  to  make  a  two-phase  flow  approach  one-phase  flow 
inside  the  entire  supersonic  nozzle  is  through  a  reduction  in  particle  mass 
fraction. 


Z/rt 


.  5.  Axial  Flow  Mach  Numbers  for  Different  Particle  Sizes 
(IUS  small  motor  exit  cone) 
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V.  NUMERICAL  RESULTS — THREE-DIMENSIONAL  NOZZLES 


The  nozzle  cross  section  considered  in  this  study  is  given  by  the 
superelliptic  equation6 


W 


for  n  >  2 


(11) 


The  variation  of  the  geometric  parameters  along  the  longitudinal  z  direction 
is  described  by  a  cosine  curve  of  expansion  and  compression  profile  of  the 
following  form: 


he  +  hi  **  “  h^  /z  -  z  \ 

h  *  - ; -  +  - = -  cos  ( -  it  J  for  z.  <  z  <  z 

2  2  \  za  /  1  "  ~  8 

*  h  for  z  >  z  ( 12 ) 

e  —  s 

where  h  stands  for  r  ,r.,  or  n,  and  z  stands  for  some  location  between 

a’  b’  s 

the  initial  z^  and  the  end  zg  station  along  the  longitudinal  z  direction. 

In  this  study,  z.  -0.0  and  z  ■  5.0.  The  reference  scale  L  in  the  two- 

X  3 

phase  flow  calculation  has  been  set  equal  to  unit  foot,  and  the  initial  cross 
section  at  z  ■  0.0  is  a  circular  profile  with  v ^  ■  rb  ■  1  and  n  *  2.  The 
flow  at  the  inlet  plane  is  a  uniform  Mach  3. 

Figure  6  illustrates  the  cross-sectional  grids  (21  x  19,  21  points  in 
r  direction  and  19  in  6  direction)  at  different  z  stations  for  a  rounded 
square  nozzle  with  a  90°  cross-sectional  flow  region.  The  superelliptic 
parameters  vary  from  the  initial  values  to  rg  ■  rfa  ■  1.5  and  n  ■  5  at 
zg  *  5.0,  according  to  Eq.  (12).  Figure  7  shows  the  computed  three- 
dimensional  Mach  number  contours  for  a  clean  gas  flow  (y  *  1.4),  where  the 
index  N  is  set  to  1  in  Eq.  (1). 


23 


Fig.  6.  Cross-Sectional  Grids  for  Rounded  Square  Nozzle 


Fig.  7.  Mach  Number  Contour  for  Rounded  Square  Nozzle  (one-phase  flow) 


For  two-phase  calculation,  the  gas  and  particle  data  are  taken  to  be  the 
same  as  those  used  in  Ref.  3. 

Gas  Phase  Particle  Phase 

5  -  2.68  KJ/kg-K  (0.64  Btu/lb  -°R)  C.  -  1.38  kJ/kg-K  (0.33  Btu/lb  -°R) 

p  m  j  m 

M..  -  8.88xl0"5Pa.s  (5.97xl0“5lb  /ft-s)  m.  -  3203.69  kg/m3  (200  lb  /ft3) 
ci  oi  j  m 

Pr  -  0.45  0  -  28. 8Z 

A  -  0.664 

The  two-phase  index  N  is  set  to  2,  and  unit  velocity  lag  and  temperature 
ratio,  A  »  XT  ■  1,  is  assumed  at  the  inlet  plane,  z  ■  0.  The 
computed  Mach  number  contours  for  a  two-phase  flow  (y  ■  1.4,  Fj  *  6  u) 
are  given  in  Fig.  8. 

Figure  9  compares  the  boundary  Mach  number  distribution  for  various  one- 
and  two-phase  flows  along  the  meridional  direction  at  z  *  2,  4,  6,  and  8,  and 
Fig.  10  presents  the  similar  results  along  one  of  the  symmetry  boundary 
planes.  Because  of  flow  symmetry  with  respect  to  a  0  *  45°  plane  for  the 
rounded  square  nozzle,  only  a  45°  region  is  plotted  in  Fig.  9.  At  the  same 
gas  specific  heat  ratio,  the  one-phase  flow  has  higher  Mach  numbers  on  the 
boundary  than  that  of  the  two-phase  flow,  except  near  the  region  downstream  vi 
the  shock  as  evidenced  at  z  ■  8  in  Pig*  10,  since  the  imbedded  shock  strength 
is  stronger  for  the  clean  gas  flow  than  for  the  flow  with  particles.  The  wave 
structure  in  the  two-phase  flow  is  different  from  that  of  the  one-phase  flow. 
For  both  the  one-  and  the  two-phase  flow  with  the  same  particle  size,  lower 
gas  specific  heat  ratio  results  in  lower  Mach  numbers.  Since  hot  combustion 
product  from  a  solid  rocket  chamber  usually  contains  gaseous  phase  with  lower 
specific  heat  ratio  than  that  of  the  gas  in  the  cold  flow  test  facility,  the 
data  obtained  from  the  cold  flow  test  should  be  used  with  discretion  when 
utilized  for  the  flight  design  consideration.  Moreover,  at  the  same  gas 
specific  heat  ratio  and  particle  mass  fraction,  the  gas  phase  in  large  particle 
two-phase  supersonic  flow  is  more  susceptible  to  the  nozzle  boundary  geometry 
change  Chan  Chat  in  small  particle  flow,  implying  that  the  small-sized  particle 
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ONE-PHASE,  y  =  1.4 


-  TWO-PHASE,  y  =  1.4/7j  =  6.0^ 

- TWO-PHASE,  y  =  1.4,77  =  1.0  M 

.  TWO-PHASE,  y  =  1.19,77  =  6.0  ^ 

- ONE-PHASE,  y  =  1.19 


8  (deg) 


Fig.  9.  Boundary  Mach  Number  Distribution  Along  Meridional  Direction 
for  Rounded  Square  Nozzles 
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Fig.  10.  Boundary  Mach  Number  Distribution  Along  x  Direction 
for  Rounded  Square  Nozzle 
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acts  more  effectively  to  retard  both  gas-phase  expansion  and  compression  than 
that  of  large-sized  particles.  Figure  11  is  the  computed  three-dimensional 
particle  density  contour  at  y  1  1.4,  r,  *  6  A  clear  particle-free 
zone  is  observed  from  the  calculated  results.  The  particle  impinges  on  the 
boundary  at  z  a  6.5,  indicating  a  region  of  high  heating  and  erosion  inside 
the  nozzle.  The  information  is  useful  to  nozzle  insulation  study. 

Whether  a  relaxation  zone^  appears  behind  the  imbedded  shock  wave  in 
a  nonuniform  multidimensional  two-phase  nozzle  flowfield  depends  upon  the  flow 
condition  in  front  of  the  shock  wave  and  upon  the  particle  size.  The  relaxa¬ 
tion  zone  is  regarded  to  be  the  transition  region  where  particle  velocity 
immediately  behind  the  shock  is  higher  than  gas  velocity,  and  particle  phase 
reaches  equilibrium  state  with  gas  phase  further  downstream  of  the  shock. 

Figure  12  shows  the  gas  and  particle  velocity  distribution  along  one  of  the 
synsnetry  boundary  planes  at  z  m  6  and  z  ■  8  for  the  flow  with  y  *  1.4,  r^ 

■  lp.  For  this  small  particle  two-phase  flow,  the  particle  velocity  downstream 
of  the  shock  is  higher  than  gas  velocity  and  a  relaxation  zone  ensues,  similar 
to  that  discussed  in  Refs.  14  and  15.  However,  for  a  large  particle  flow, 
different  shocked  flow  behavior  is  observed.  Figure  13  shows  the  gas  and 
particle  velocity  distribution  for  the  flow  with  r .  *  6  |j  at  the  same 
locations.  No  relaxation  zone  is  possible  for  this  large  particle  flow,  since 
both  the  particle  velocities  in  front  of  and  behind  the  imbedded  shock  wave 
are  lower  than  the  corresponding  gas  velocities  and  an  equilibrium  state  for 
gas  and  particle  does  not  exist.  Obviously,  the  results  shown  in  Figs.  12  and 
13  are  the  direct  consequence  of  assumed  particle  drag  coefficient  and  heat 
transfer  formulation.  Different  formulations  would  affect  the  computed  flow- 
field.  The  discussions  given  in  Refs.  14  and  15  on  the  existence  of  a 
relaxation  zone  behind  a  shock  wave  in  a  two-phase  flow  are,  strictly 
speaking,  applicable  only  for  a  homogeneous  uniform  one-dimensional  two-phase 
mixture  in  front  of  the  shock  wave. 
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Fig.  11.  Particle  Density  Contour  for  Rounded  Square  Nozzle 


X 


Fig.  12.  Gas  and  Particle  Velocity  for  r:  »  ly — with 
Relaxation  Zone  Behind  Imbedded JShock  Wave 
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Fig.  13. 


Gas  and  Particle  Velocity  for  rj  ■  6p — without 
Relaxation  Zone  Behind  Imbedded  Shock  Wave 
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Figure  14  illustrates  computed  flow  Mach  number  at  the  centerline  of 
the  symmetry  plane.  The  similar  effect  of  particle  size  variation  on  the  flow 
behavior  as  that  of  axisymmetric  flow  shown  in  Fig.  5  is  observed  for  the 
three-dimensional  flow  inside  the  rounded  square  nozzle.  The  computer  program 
with  the  array  dimensioned  for  a  31  x  37  cross-sectional  grid  occupies  151,000 
(octal)  small  core  and  214,000  (octal)  large  core  on  a  CDC  Cyber  176  computer, 
and  a  typical  run  for  the  two-phase  flow  inside  the  rounded  square  nozzle  with 
21  x  19  cross-sectional  grid  takes  2  min,  10  sec. 

The  cross-sectional  grids  at  different  z  stations  inside  an  elliptic 
nozzle  with  90°  cross-sectional  flow  region  are  illustrated  in  Fig.  15.  The 
geometric  parameters  vary  smoothly  from  the  initial  values  to  r^  ■  1.0,  r^ 

*  1.5,  and  n  *  2  at  z  *  5.0.  The  same  inflow  condition  as  the  rounded  square 
nozzle  is  used  for  the  elliptic  nozzle  calculation.  Figure  16  depicts  the 
computed  particle  density  contour,  which  shows  that  a  smaller  particle-free 
zone  existed  in  the  elliptic  nozzle  than  that  in  the  rounded  square  nozzle. 
Figures  17  and  18  compare  the  difference  in  Mach  number  contours  in  the 
elliptic  nozzle  for  one-  and  two-phase  supersonic  flows.  The  compression 
waves  generated  from  the  boundary  coalesce  into  a  weak  shock  wave,  which 
impinges  on  the  centerline  of  the  elliptic  nozzle  earlier  than  that  of  the 
rounded  square  nozzle  in  both  one-  and  two-phase  flows.  The  comments  made  in 
Ref.  6  regarding  the  dependence  of  the  flowfield  on  the  nozzle  geometric 
configuration  for  the  clean  gas  flow  are  equally  applicable  to  the  two-phase 
flow. 

Finally,  a  hybrid  nozzle,  which  is  composed  of  half  circular  (0  ^  0  £  90°) 
and  half  rounded-square  (90°  £  0  <_  180°)  regions  is  depicted  in  Fig.  19  to 
illustrate  the  calculation  with  a  180s  cross-sectional  flow  region  (21  x  25). 
The  auperelliptic  parameters  vary  from  the  initial  values  Co  ra  ■  r^  ■  1.5 
and  n  ■  2  for  the  circular  portion  and  to  r  ■  r.  *  1.5  and  n  ■  5  for  the 
rounded-square  portion  of  Che  nozzle  at  z  1  5.0.  Figure  20  shows  the  computed 
particle  density  contours  for  the  two-phase  flow  with  r ^  *  6y  and  y  m  1.19. 
Because  of  a  smaller  gas  specific  heat  ratio,  Che  value  of  Che  nondimens ional- 
ized  particle  density  given  in  Fig.  20  is  smaller  than  that  in  Fig.  11.  A 


Fig.  16.  Particle 


Fig.  17.  Mach  Number  Contour  for  Elliptic  Nozzle 
(two- phase  flow) 


Fig.  IB.  Mach  Number  Contour  for  F.lliptic  Nozzle 
(one-phase  flow) 


Cross-Sectional  Grids  for  Hybrid  Nozzle 


Fig.  20.  Particle  Density  Contour  for  Hybrid  Nozzle 


uniform  Mach  3  flow  at  Che  inlec  plane  means  chat  changing  y  from  1.4  to 
1>19  causes  Che  reduction  of  Che  ratio  of  gas  static  to  reference  stagnation 
density  by  a  factor  of  1.97,  according  to  the  isentropic  relationship.  The 
same  factor  is  carried  over  to  the  dimensionless  particle  density  value  of  the 
inlet  plane,  since  the  particle  mass  fraction  and  the  reference  state  remain 
the  same.  Figures  21  and  22  compare  the  Mach  number  contours  for  the  one-  and 
two-phase  flows,  respectively.  It  is  obvious  from  these  Mach  number  contour 
curves  that  the  introduction  of  rj  ■  6y  particles  in  the  gas  flowfields 
alleviates  the  contour  clustering  and,  hence,  the  shock  strength  by  approxi¬ 
mately  30  percent.  This  perhaps  explains  why  the  performance  calculation  from 
the  method  of  characteristics  fortuitously  predicts  a  result  which  agrees 
relatively  well  with  that  from  motor  test  firing,  even  though  a  shockless 
flowfield  is  tacitly  assumed  to  prevail  in  Ref.  4. 


Fig.  21.  Mach  Number  Contour  for  Hybrid  Nozzle  (two-phase  flow) 


Fig.  22.  Mach  Number  Contour  for  Hybrid  Nozzle  (one-phase  flow) 


VI.  CONCLUSION 


Idealized  two-phase  supersonic  flows  inside  an  axisymmetric  IUS  small 
motor  nozzle  and  various  three-dimensional  nozzles  have  been  studied 
numerically.  The  study  stressed  the  importance  of  the  proper  treatment  of 
boundary  flow  variables  for  a  two-phase  flowfield  calculation,  even  in  the 
absence  of  a  pertinent  particle  impingement  model.  Subject  to  the  provision 
that  the  particular  particle  drag  law  for  spherical  particles  and  the  simple 
particle  convective  and  radiative  heat  transfer  relationships  have  been 
adopted  in  the  study,  the  computed  results  revealed  the  three-dimensional 
two-phaae  supersonic  nozzle  flow  structures,  which  differed  noticeably  from 
those  of  the  clean  gas  flows.  The  study  also  emphasized  the  existence  and 
complexity  of  the  imbedded  shock  inside  the  three-dimensional  supersonic 
nozzles.  Due  to  the  restriction  of  practical  nozzle  physical  dimensions,  the 
complex  imbedded  shock  structure,  often  ignored  or  assumed  nonexistent  in  the 
calculation  with  some  simplified  analysis,  is  almost  ubiquitous  in  a 
three-dimensional  supersonic  nozzle  and  can  be  computed  with  the  present 
technique.  Obviously,  a  three-dimensional  one-  or  two-phase  supersonic  nozzle 
flow  generally  would  have  nonuniform  flow  at  the  inlet  plane,  and  the  starting 
data  should  be  obtained  from  a  three-dimensional  subsonic-transonic  flow 
solution,  which  is  the  subject  of  follow-on  studies. 
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n  l  n  I 


NOMENCLATURE 


a. ai 

A 

Ai 

b, b  i 

BJ 

°D 


e 

e.f.g.h 

E.F.G.H 


8c**r 


h . hg , 


M 

“j 

n 


N 


Nu- 


P 

Pr 


constants  used  in  Eq.  (8) 

exponent  in  the  viscosity-temperature  Eq.  (7) 
dimensionless  friction  term  defined  in  Eq.  (2) 
constants  used  in  Eq.  (8) 

dimensionless  energy  exchange  term  defined  in  Eq.  (3) 

particle  drag  coefficient 

particle  heat  capacity 

gas  specific  heat  at  constant  pressure 

dimensionless  gas  total  energy/unit  volume 

vectors  defined  in  Eq.  (1) 

vectors  defined  in  Eq.  (9) 

momentum  transfer  parameter  defined  in  Eq.  (4) 
convective  and  radiative  parameters  defined  in  Eq.  (3) 
parameters  used  in  Eq.  (12) 

dimensionless  particle  total  energy/unit  volume 
constant  used  in  Eq.  (8) 

reference  length  scale,  e.g.,  unit  foot  or  unit  meter 
gas-phase  Mach  number 
particle  mass  density 

superelliptic  geometric  index  defined  in  Eq.  (11) 

index:  N  ■  1,  one-phase  gas-only  flow;  N  ■  2,  two-phase  flow 

particle  Nusselt  number 

dimensionless  pressure 

gas-phase  Prandtl  number 


w  -  •  •  **r 
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NOMENCLATURE  (Continued) 


q.qj 


ra*  rb 


rj 

rW 


rc 


Re; 


T,T 


j 


u,v,w 


vraaxl 

x.y.z 

Y 

6 

£e 

C 

n 

0 

*q 

M 

% 

»*tl 


dimensionless  gas  and  particle  speed,  respectively 
dimensionless  radial  coordinate 

dimensionless  nozzle  wall  radii  on  x,y  axes,  respectively, 
used  in  Eq.  (11) 

particle  radius 

dimensionless  nozzle  wall  radial  coordinate 

dimensionless  nozzle  throat  radial  coordinate 

particle  Reynolds  number  defined  in  Eq.  (6) 

dimensionless  gas  and  particle  temperature,  respectively 

dimensionless  gas-phase  velocity  component  along  z,r,6 
directions,  respectively 

dimensionless  particle-phase  velocity  component  along  z,r,@ 
directions,  respectively 

adiabatic  maximum  speed  evaluated  at  the  inlet  plane 
dimensionless  Carcesian  coordinates  shown  in  Fig,  1 
gas  specific  heat  ratio 

geometric  index;  6  *  0,  rectangular  coordinates;  6  ■  1, 
cylindrical  coordinates 

effective  emissivity  between  gas  and  particle 
transformed  dimensionless  coordinate  in  radial  direction 
normalized  radial  coordinate  ■  r/r^ 

meridional  angle  in  cylindrical  coordinates,  see  Fig.  1 

initial  velocity  lag,  ■  qj/q 

initial  temperature  ratio,  Af  ■  T/Tj 

micron,  10”^m 

gas  viscosity 

gas  viscosity  at  stagnation  state 


P 
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NOMENCLATURE  (Concluded) 


£  transformed  dimensionless  coordinate  in  meridional  angular 

direction 

p,Pj  dimensionless  gas  and  particle  density,  respectively 

0  Stefan-Boltzmann  constant 

T  (Y-l)/(2y) 


e 

(i) 


particle  mass  fraction 


Cj/Op 


Superscript 


(-) 

(-) 

Subscript 

e 

8 

i 

j 

maxi 

tl 


vectored  quantity 
dimensioned  quantity 

exit  station 
gas  phase 

initial  station  at  inlet  plane 
particle  phase 

maximum  state  evaluated  at  the  inlet  plane 
stagnation  state  evaluated  at  the  inlet  plane 


